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A complete turbulence model, where the turbulent Prandtl and Schmidt numbers are 
calculated as part of the solution and where averages involving chemical source terms are 
modeled, is presented. The ability of avoiding the use of assumed or evolution Probability 
Distribution Functions (PDF’s) results in a highly efficient algorithm for reacting flows. The 
predictions of the model are compared with two sets of experiments involving supersonic 
mixing and one involving supersonic combustion. The results demonstrate the need for 
consideration of turbulence/chemistry interactions in supersonic combustion. In general, 
good agreement with experiment is indicated. 

I. Introduction 

Accurate prediction of flows in scramjet engines requires the development of turbulence models that 
calculate the turbulent Prandtl, Pr t , and Schmidt, Set, numbers as part of the solution, and account for 
turbulence/chemistry interactions. Traditional turbulence models that only address velocity fluctuations 
have no mechanism for incorporating turbulence/chemistry interaction and require the specification of both 
Pr t and Sc t . Such numbers have a profound influence on flow predictions: a low value of Sc t can result in 
engine unstart, while a higher value may result in flame blow-out. 1 * On the other hand, Prt has an important 
effect on mixing at high speed flows. It is shown in Ref. 2, which considered the role of variable turbulent 
Schmidt number on the mixing of supersonic streams, that a value of Pr t = 0.9 gave the best fit for data 
from the experiment of Cutler et al, 3 while a value of 0.5 gave the best fit for the experiment of Burrows 
and Kurkov. 4 

In an attempt to address this problem, a series of step-by-step investigations were carried out to develop 
a model that calculates Prt and Set as part of the solution and addresses turbulence/chemistry interactions. 
Thus, in Ref. 2 the role of variable Set on supersonic mixing was considered, while in Ref. 5 the role of 
variable Prt on heat flux in the presence of shock wave/boundary interactions was examined. In a more 
recent investigation, 6 the variable Sc t formulation of Ref. 2 was extended to address reacting flows while 
assuming a fixed Pr t . 

The turbulence/chemistry interaction in Ref. 6 was studied using the multi- variate /3-PDF for mass frac- 
tions developed by Girimaji. 7 A comparison of assumed and evolution PDF’s in flows involving supersonic 
combustion by Baurle et al 8 showed that both formulation yielded comparable mean flow predictions. How- 
ever, assumed PDF’s were unable to predicted higher order correlations, such as terms involving chemical 
production source terms, with any reasonable accuracy. Similar results were encountered in Ref. 6. It is 
shown there that the use of Girimaji ’s PDF has a highly dissipative effect on the concentration variance 
resulting in poor agreement with experiment. Computations employing evolution PDF’s are time consuming 
and require excessive storage. Because of this, all terms involving production terms are modeled in this 
work. 

* Research Assistant Professor, Mechanical and Aerospacing Engineering, Member AIAA 
t Professor, Mechanical and Aerospace Engineering, Fellow AIAA 

J Aerospace Engineer, Hypersonic Airbreathing Propulsion Branch, Senior Member AIAA 


1 of 13 


American Institute of Aeronautics and Astronautics Paper 2006-0128 



The model is used to predict the flows in two sets of mixing experiments, 3,4 and the reacting experiment 
of Ref. 4. In general, good agreement is indicated. 

II. Formulation of Problem 


A. Governing Equations 

A variable Pr t and Set formulation requires equations for the variance of enthalpy and its dissipation 
rate, and the variance of concentrations and its dissipation rate. These equations were derived in Ref. 2 
and Ref. 5 for non-reacting flows. The formulation of Ref. 2 was extended in Ref. 6 to reacting flows while 
keeping Prt constant. 

The approach that has been used to derive the final set of equations for variable Pr t and Sc t in the 
presence of reactions follows the same procedure used in Refs. 2, 5, 6 and 9. This entails deriving the exact 
equations that govern the variances of concentrations and enthalpy and their dissipation rates from the 
Navier-Stokes equations and model the resulting equations term by term. This insures that relevant physics 
is incorporated into the model. Dimensional and tensorial consistency, Galilean invariance, coordinate system 
independence, and absence of wall or damping function characterize the resulting set of equations which are 
given in the Appendix. 

B. Turbulence/Chemistry Interaction 

The equation for mass fraction variance, ay, contains the term: 

where is the fluctuation of the mass fraction of species to, and uj m is its mass production rate. Similarly, 
the equation that governs the enthalpy variance contains the term 

h"d) m A/ty m 

where h" is the enthalpy fluctuation and A hf m is the heat of formation of species m. 

Traditionally, the above terms are evaluated by using an assumed or evolution PDF’s or ignored com- 
pletely. The assumed PDF’s are usually a product of Girimaji’s multi-variate /3-PDF for mass fraction 
fluctuations and a Maxwellian for temperature fluctuations. Comparisons of the predictions of assumed and 
evolution PDF’s have been conducted by Baurle et al 8 on supersonic combustion of parallel stream. It was 
shown there that both formulations give comparable mean values. However, assumed PDF’s were unable to 
produce the correct values of the higher order correlations. In particular, they gave the wrong sign for cor- 
relations involving mass production rates. This is why, in the absence of evolution PDF’s, better predictions 
are obtained when correlations involving mass production rates are set to zero. 

As will be shown, in the Results and Discussion section below, setting correlations involving mass produc- 
tion rate to zero is not an option for the current formulation. Because evolution PDF’s require an excessive 
amount of time and storage, these terms are modeled here. Thus, 

2Y / Y^JZ = CysY / \/W 2 ^ ( 1 ) 

m m 

and 

E^ m Ahf t m — Ch,12 \! h" 2 ^ ] dj m Ahf tm (2) 


C. Numerical Procedure 

A modification of REACTMB, 10 a code that has been under development at North Carolina State 
University over the last several years, is employed in this investigation. It is a general purpose parallel 
Navier-Stokes solver for multi-component multi-phase reactive flows at all speeds. It employs a second 
order essentially non-oscillatory (ENO) upwind method based on Low Diffusion Flux Splitting Scheme of 
Edwards 11 to discretize the inviscid fluxes while central differences are employed for the viscous and diffusion 
terms. Planar relaxation is employed and the code is parallelized using domain decomposition and message 
passing(MPI) strategies. 
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D. Model Constants 


The model constants developed in Refs. 2, 5 and 6 remain unchanged. The final set of model constants 
are summarized in Table 1 for concentration variance and its dissipation rate, and in Table 2 for enthalpy 
variance and its dissipation rate. 

Table 1. Model constants for ay and ey equations 


C Y 

C'y.i 

Cy, 2 

Cy, 3 

Cy, 4i 

Cy, 42 

Cy, 5 

Cy, 6 

Cy, 7 

Cy, 8 

Cy, P 

Cy, 9 

CTh 

0.065 

1.0 

0.095 

-0.025 

0.45 

-1.0 

1.0 

0.5 

0.78125 

0.25 

-0.1 

1.0 

0.5 


Table 2. Model constants for h" 2 and equations 


c h 

Ch, 2 

c hA 

C h , 5 

C h , 6 

Ch, 7 

C h ,s 

C h , 9 

Ch, 10 

C h , n 

Ch, 12 

0.0648 

0.5 

-0.4 

-0.04 

-0.12 

1.45 

0.7597 

0.87 

0.25 

-1.5 

-0.75 


III. Results and Discussion 

A theory that is developed to predict Pr t and Sc t as part of the solution should apply for both reacting 
and non-reacting flows. Because of this, the present theory is validated by two sets of experiments involving 
supersonic mixing, 3,4 and one experiment involving supersonic combustion. 4 In the experiment of Cutler 
et al, a coaxial nozzle was designed to produce two uniform coaxial jets at exit. The center jet consists 
of 95% of He and 5% O 2 by volume at a Mach number M = 1.8, while the outer jet is air at M = 1.8. 
A schematic of the experiment is shown in Fig. 1. The grids employed are the ones used in Refs. 2 and 3. 
The fine grid consists of 188, 080 cells and is decomposed into 13 blocks for parallel computing while the 
intermediate grid deletes every other point in the axial direction. All results employed here model the flow 
in the nozzle, employ the fine grid and use the axisymmetric version of REACTMB. 

The second set of experiments are those of Burrows and Kurkov. 4 A schematic of the experiment is shown 
in Fig. 2. Hydrogen is injected into the test section through a nickel injector parallel to the vitiated main 
flow. The mixing case employed nitrogen in place of air for the main flow. At the entrance of the test section, 
M = 2.44, the static pressure is one atmosphere, and the static temperature is in the range 1250-1270 K for 
the reacting case, and about 1150 K for the mixing case. In both cases, hydrogen was injected at M = 1, 
matched pressure and a total temperature slightly above the ambient temperature. The two grids that are 
employed here are those used in Ref. 6. Each grid consists of 15 blocks. The first grid has 86,643 cells, while 
the second has 104,428 cells. The fine grid reflects grid refinement in the blocks where mixing of the two 
streams takes place. Rather than using measured conditions at the inlet of the test section, the flows in both 
hydrogen and nitrogen/air nozzle were computed. It was necessary to iterate on inflow conditions of both 
nozzles to arrive at the stipulated conditions at the exit of each nozzle. All results presented here employed 
the fine grid. 

In Ref. 2, calculations of Cutler et al experiments 3 employed a variable Sc t , and a Pr t = 0.9. Results are 
presented in Figs. 3-5, which compare predictions of current theory at selected stations with those of Ref.2 
and the experiment. As is seen from Fig. 3, the results of Ref. 2 for mass fraction of He at x = 261mm are 
better than the current prediction. Figures 4 and 5 show that predictions for velocity and Pitot pressure are 
comparable. 

Figure 6 shows a comparison of mass fraction prediction with the mixing experiment of Ref. 4. As is seen 
from the figure, good agreement is indicated. Similar results were obtained in Ref. 6. 

Two sets of figures are presented for the reacting case of Ref. 4. For this calculation, the seven-species, 
seven-reaction H 2 -AI 1 : model used in Ref. 6 which was originally developed by Jachimowski 12 is employed. In 
the first, terms involving averages of chemical source terms are ignored, while in the second, the contributions 
of these terms are included. As is seen from Fig. 7, poor agreement with the experiment is indicated. Figures 8 
and 9 show contours of Set and Prt . It appears that the main cause of the discrepancy is a result of a reduced 
Prt near the mixing region. This has the tendency of promoting heat transfer and early combustion. Figure 
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10 shows that when the contributions of the chemical source terms are included, much better agreement 
with the experiment is indicated. Figures 11 and 12 show contours of Set and Prt, respectively. 

Based on the above, two relevant observations can be made. The first is that modeling of averages of 
terms involving chemical source terms is a viable option. When this approach is compared with approaches 
requiring assumed or evolution PDF’s, a great deal of computational efficiency is achieved. Second, a 
relatively inexpensive calculation of variable Prt and Set can be obtained by assuming a value for the Lewis 
number, and eliminating either the equations for enthalpy variance and its dissipation rate, or those for the 
variance of concentrations and its dissipation rate. This, however, will result in ignoring one of the averages 
involving chemical source terms. Because inclusion of such terms is important, assuming a constant Lewis 
number is not recommended. 


IV. Conclusions 

A complete turbulence model where both Pr t and Set are calculated as part of the solution, and where 
averages of terms involving chemical source terms are modeled is presented. Thus, calculations of turbulent 
flow become similar to those of laminar flows in the sense that all that is required is to specify initial and 
boundary conditions. 

Because the resulting algorithm does not require the use of an assumed or evolution PDF, it is compu- 
tationally efficient, especially when one deals with complex three-dimensional geometries characteristics of 
proposed scramjet designs. 

The resulting algorithm, which is based on the exact Navier-Stokes equations, is dimensionally and 
tensorially consistent, Galilean invariant, coordinate system independent and free of damping and wall 
functions. Although, the algorithm is applied to relatively simple geometries, past experiences suggest that 
such an approach works well for complicated geometry without having to adjust any of the model constants. 

Finally, averages involving chemical source terms are important and shall always be included in combus- 
tion calculations. 
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A. Model Equations 


The Favre-averaged species conservation equations can be written as 


d _ n dY m 


where 


;{pYm) + g x _ (pUjY m ) - dX ' I pD dx _ pYyi 1 + w m . 


~p y ™ u " = P Dt yr:' 


p is the density, Y m is the mass fraction of species m, Ui is the velocity, D is the laminar diffusion coefficient 
and D t is the turbulent diffusion coefficient. 

The variance of mass fractions, cry, is defined as 




while its dissipation rate, ey, is defined as 




The cry-equation is given as 


where 


d 9 d I" day 

g,(p-y) + ^ [«(D + cy,D,)— 


+2 ^ pD t — — j - 2 pey + 2 ^ LOmXrn 


+ C Y *jy ma.^ — ,0.0 


n 1 (n k °Y v t 1 k 2 

Dt = ^(CV — + — )» > Vt = C nY7 

2 ey a h 2 

2 J2XnK = Cy,sY,\f^ 2Z ^’ 


P is the pressure, i p is the turbulent kinematic viscosity, k is the turbulent kinetic energy and £ is the 
enstrophy. 

The ey- equation is given as 


wX Y)+ -Ly v) = 


(l din duA _ <9 yr 2 dY m 

+2pev ( + CY ’ 2bjk XA + c Y3P k 22xrv Y m XX 


T pDCy/iiDf 'yA 


dxjdxj 


-n^Y , 42 . / V "2 ^ 

■pD XV Y ™ a — T~ 

Ty ^ v OXkOXk 


- n Cy,6 I ^ I y. I Y^ , /w" 

pDt V y — — Cyjp h X\ Y m 

Ty z -—' \ OXj / Ty Ty z -—' v 

m. \ J / m 
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where 


7 Tjk | 2 

bjk ~^k + 3 6jk ’ 


Tjk — P u j u k 


The turbulent Schmidt number is defined as 



(A.7) 


The mean energy equation can be written as 


| ( pi) + = w “ t; + * “ ^ ( p ' 1 " 1 


dt 


Dt dx; 


where 


ft = 




— ph"u '■ = 


tin — 


— -/ ^ i n 7 'i 

ft,i-PK a +AL^m g ) 

J m J 

r dUi 

t »ai] + pt ’ C = K 


O /c _ 1 dUi 0Uj. 

3% aa , fc )- 2 ( a Cj - + dx? 


h is the enthalpy, qt is the laminar heat flux, and a t is the turbulent diffusivity. 
The enthalpy variance ( h " 2 ) equation can be written as 


i(p h "V2) + ^(pu j h"y2) = ±- 


dt 


dxj 


p{"1oc + a t C h , 2 ) 


dh" 2 /2 

dXn 


+2/ryS». 


d ,q t i s d ,q* is 


[ dxj p dxi p 

P h " 2 JT 

OX 

- s ^h"Cj m A hf. 


d L/X j fJ 


cii) ■ f) h 

-(7 - l ^P h " 2 ~Q^. ~~ qt ’ i faP- + 2C M7 P\J ^" 2 C - lP e h 


with 


7 

= C p /c v , 


a t 

= 2 ( C'hkxh + 

0.89/ 7 

Th 

)§ f-e 

(-£S 1 ^ 

II 



/ dh " \ 2 
“ l i>x, J ’ 



= C'fc.iaV^] 



h"<j m Ahf >r 

m m 

where is the dissipation rate of the enthalpy variance, and a is the laminar diffusivity. 


(A.8) 


(A.9) 
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The equation for the dissipation rate of enthalpy variance is taken as 


d 


d 


= -pe h ( C h , 5 b jk - — ) ^7 


dt 


Ojk 


du.j 




dVhP 2 


dh 


dxj dxj 


dx k 

d 

dx; 


(ja + C h pat) 


de h 

dxj 


c< 

^h,8 

+Ch,ll^h 


Qt.j dh - ( Ch 9 

- ipe-h 1 


C h . 


10 


T h OXj 

Dp 
Dt 


\ Xfj T k 

(DP 


P 

-=■ max 
P 


v Dt 


,0.0 


where 


Tk = 




The model constants, Ch and Ch, 1 - 12 , are given in Table 2. The turbulent Pr* is defined as 

T3 V t 

Pr t = — 
a t 



Figure 1. 


Schematic of Experiment Setup (Ref. 3) 
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Figure 2. Schematic of Experiment Setup (Ref. 4) 



Figure 3. Comparison of computed and measured He-O 2 mass fraction(Ref. 3) 
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Figure 4. Comparison of computed and measured velocities (Ref. 3) 



Figure 5. Comparison of computed and measured Pitot pressure (Ref. 3) 
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Figure 6. Comparison of computed and measured volume fractions (Ref. 4), mixing case 



Figure 7. Comparison of computed and measured volume fractions (Ref. 4), reacting case, without chemical 
source terms 
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Figure 9. Prandtl number contours, reacting case, without chemical source terms 
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Figure 10. Comparison of computed and measured volume fractions (Ref. 4), reacting case, with chemical 
source terms 
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Figure 11. Schmidt number contours, reacting case, with chemical source term 
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Figure 12. Prandtl number contours, reacting case, with chemical source term 
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